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ABSTRACT 


An  orbit  averaging  numerical  integration  method  for 
determining  the  variation  in  shape  and  orientation  of  an  earth 
satellite  orbit  is  discussed. 
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ORBIT  AVERAGED  VARIATION  OF 
EARTH  SATELLITE  ORBITAL  ELEMENTS 


1.  INTRODUCTION 

In  this  note  we  sketch  the  derivation  of  the  equations  for  the  osculating 
elliptic  orbital  elements  of  an  earth  satellite  in  terms  of  an  angular  indepen¬ 
dent  variable  and  discuss  an  orbit  averaging  technique  of  integrating  the 
equations.  We  employ  equations  which  are  valid  for  small  eccentricity  and 
inclination. 

Because  we  are  interested  in  earth  satellite  orbits  for  which  moon 
perturbations  are  difficult  to  handle  analytically,  we  discuss  the  numerical 
integration  of  the  orbit  averaged  equations. 

Gaussian  quadrature  of  orbit  averaged  definite  integrals  is  much  faster 
in  computer  time  than  the  exact  integration  of  the  equations  of  motion.  Closed 
form  analytic  formulas  for  the  definite  integrals  would  be  even  faster  to  evalu¬ 
ate,  but  as  was  said  the  derivation  of  such  formulas  is  extremely  intricate 
for  far  out  satellites. 

The  method  of  orbit  averaging  yields  reasonably  accurate  predictions 
of  the  shape  and  orientation  of  the  orbit,  but  not  of  the  specific  position  in 
the  orbit.  A  first  order  method  of  evaluating  the  orbit  averaged  definite  inte¬ 
grals  might  cause  the  predictions  to  get  out  of  phase  with  time.  A  second 
order  method  is  described  which  uses  a  first  order  orbit  averaged  linear 
variation  in  the  elements  in  evaluating  the  definite  integrals.  It  would  be 
more  accurate  to  use  a  first  order  analytic  non-orbit  averaged  variation  in 
the  elements  in  the  second  order  evaluation  of  the  definite  integrals.  Any 
such  formulas  would  be  easy  to  add  to  the  computer  program  that  is  being 
written  to  implement  the  formulas  in  this  note. 

Explicit  formulas  for  the  perturbing  accelerations  due  to  the  sun,  moon 
and  earth  gravitational  potential  harmonics  are  presented  in  order  to  docu¬ 
ment  the  computer  program  as  completely  as  possible. 
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II. 


CARTESIAN  EQUATIONS  OF  MOTION 


The  equations  of  motion  of  an  earth  satellite  are 


where 


^2  k 

d  X 


dt 


k 

X 


k 

X 


-  !^  +  f’' 

r 

k  dx^  k+3  , 

X  ,  =  X  when  t 

o  dt  o 


2,  3  (1) 


k  cartesian  position  component  of  satellite 
relative  to  earth  (k  =  1,2,3) 


X  =  k  cartesian  velocity  component  of 

satellite  relative  to  earth  (k  =  1,2,3) 


r 


1  TT 

Z  (  ) 

j=l 


(JL  =  gravitational  constant  times  mass  of  earth 

3  -2 

(units  are  E  T  ) 
k  th 

F  =  k  component  of  acceleration  on  satellite  additional 
to  the  central  force  acceleration  (k  =  1,2,3) 


If  the  perturbing  accelerations  (F  ,  F  ,  F  )  were  zero,  the  motion 
satisfying  (1)  would  follow  an  elliptic  trajectory.  Even  with  non- zero  perturb¬ 
ing  accelerations  there  are  elliptic  elements  associated  with  the  postion  and 
velocity  (x\  .  .  .  ,  x^)  at  time  t  ,  called  the  osculating  elements.  These  are 
the  elements  of  the  elliptic  orbit  that  the  satellite  would  follow  given  the  position 
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I 


and  velocity  (x  ,  .  .  .  ,  x  )  at  time  t  if  the  perturbing  accelerations  were 
turned  off  after  time  t  . 

The  elliptic  elements  we  shall  employ  are 

a  =  Semi-major  axis  (distance  unit  consistent  with  the  distance 
unit  in  fi) 

e  =  Eccentricity  (0  ^  e  <  1) 

I  =  Inclination  (0^  ^  I  ^  180^) 

n  =  Right  ascension  of  ascending  node  (0^  ^  0  <  360^) 

rn  =  Argument  of  perigee  (0^  ^  uj  <  3  60^) 

=  Mean  anomaly  at  time  t^  (0  ^  <  3  60^) 

See  Ref.  1  or  any  book  on  celestial  mechanics  for  the  formulas  for  the  position 
and  velocity  as  a  function  of  time  in  terms  of  the  orbital  elements  and  for  the 
orbital  elements  in  terms  of  the  position  and  velocity  at  a  given  time.  Addi¬ 
tional  notation  for  elliptic  orbits  which  we  shall  employ  is 

2 

p  =  a(l  —  e  )  =  semi-latus  rectum 

1/2  -3/2 

n  -  p.  a  =  mean  motion 

M  =  +  n(t  —  t^)  =  mean  anomaly  at  time  t 

§  =  eccentric  anomaly 

Y  =  true  anomaly 

T\  =  0)  +  Y  =  angle  measured  along  the  orbital  plane  from  the 

ascending  node  to  the  true  position  of  the  satellite 
Tjy  =  ff)  +  Q  =  longitude  of  perigee 

rj  =13  +  0 

In  order  to  study  the  behavior  of  the  osculating  elements  with  time  the 
cartesian  equations  of  motion  (1)  can  be  numerically  integrated  and  the 
osculating  elements  evaluated  at  each  tabular  point.  However,  numerical 
integration  of  (1)  takes  a  fair  amount  of  computer  time  and  becomes  pro¬ 
hibitive  when  studying  a  large  variety  of  orbits. 
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III.  EQUATIONS  FOR  OSCULATING  ELEMENTS 


Equations  of  motion  (I)  can  be  written  in  the  form 


dx 

“dT 


=  X 


k+3 


dx 


k+3 


dt 


k  k 

X  =  X 

o 


+  F 


>k  =  1,  2,  3 

(2) 


k+3  k+3 

X  =  X  when  t 

o 


=  t 


Let  (b\  ,  B^ .  6^,  B^)  =  (a,  e,  1,0,  UU,  M^)  be  the  osculating  elliptic 

orbital  elements  at  time  t  with  the  functional  relationship 


k  k,„l  ,  ,  ,  / 

X  =x{B,...,B,t),k  =  l . 6 


being  given  by  the  usual  elliptic  orbit  formulas  in  such  a  way  that 


(3) 


5x 

at 


=  X 


k+3 


bx 


k+3 


px 


at 


k  =  1.  2.  3 


(4) 


We  have 


dx 

“dF 


6 

=  E 


ax  dB” 


j=i  aB^ 


dt 


ax 

at 


k  =  1,  . . . ,  6 


(5) 
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Equations  (2),  (3),  (4)  imply 


6 

Z 


j=l 


dt 


0 


6 

E 

j=l 


Sx 


k+3 


SB- 


dt 


k  =  1,  2,  3 


(6) 


We  introduce  the  Lagrange  brackets 


For  a  given  value  of  i  (between  1  and  6)  multiply  the  first  equation  in  (6) 

^^k+3  k 

by  _  - ; —  and  the  second  equation  in  (6)  by  for  k  =  1,  Z,  3  and 

SB^  SB^ 

then  add  the  resulting  six  equations.  We  obtain 


I  [e‘ ,  b']  ^ 

J=1 


3 

Z  F 
k=:l 


k  9x 


1, . .  .  ,  6 


(8) 


The  Lagrange  brackets  [B^  ,  3,re  evaluated  and  equations  (8)  solved  for 

dBV^t  (j  =  1,  .  .  .  ,  6)  in  most  books  on  celestial  mechanics;  see,  for  example, 

Ref.  2,  pages  273-307.  We  shall  just  give  the  final  result. 

Let  be  a  unit  vector  at  the  satellite  pointing  from  the  earth  to  the 

satellite;  let  be  a  unit  vector  normal  to  0^  in  the  plane  of  the  instantaneous 

position  and  velocity  vectors  of  the  satellite  and  making  an  acute  angle  with  the 

velocity  vector;  and  let  0  =  ?  x  0  .  Denote  the  components  of  the 

^  w  r  s  ^ 
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12  3  — *  — ♦ 

acceleration  (F  ,  F  ,  F  )  in  the  T  ,  T  and  e  directions  by  R,  S  and 
'  ’  r  s  w  ^ 

W  ,  respectively.  Then  Gauss'  form  for  the  equations  for  the  osculating  ele¬ 
ments  are 


da 

dt 


[R  e  sin  Y 


+ 


de 

dt 


/l  - 

na 


[  R  sin  Y 


+  S(cos  § 


+  cos  Y )] 


dt 


rW 


cos(  uj  +  Y ) 


•  T 

sin  I  ^ 


dm 

dt 


rW 


s  in  (  m  +  Y  ) 


v/l  - 

nae 


[-  R  COS 
—  cos  I 


Y 

dt 


+ 


S(  —  +  1 )  sin  Y ] 
P  ^ 


dM 

dt 


n  + 


S 

nae 


+ 


1  -  e 
nae 


2 


cos  Y 


] 


The  equation  for  the  initial  mean  anomaly  is  easily  derived  from  the 

above  if  rather  than  M  is  regarded  as  the  primary  orbital  element. 
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IV.  CHANGE  OF  INDEPENDENT  VARIABLE 


The  eccentric  anomaly  §  ,  true  anomaly  Y  and  mean  anomaly  M  in 
an  elliptic  orbit  are  related  by 


* 

tan  ^ 


/ 


1  +  e  ^  § 

- -  tan  - 

1  —  e  2 


M  =  5  ~  e  sin  ^ 

From  these  formulas  and  equations  (9)  we  find  after  some  very  intricate 
elliptic  orbit  computations 


^  ^  ^  ^  (pp)^^^  ( .  _  £_ 

dt  dt  dt  2  y  pp 


cot  I  sin  r|  •  W 


) 


For  the  compound  angular  quantity  r|  =  T)  +  we  then  have 


dn  _  (pp) 

dt  2 

r 


1/2 


r  sin  r|  »  W  1  —  cos  I 

]  .1/2  sin  I 

(J^P) 


We  define 


K  = 


1  +  — 


3 


(10) 


1  —  cos  I  . 

- : - z -  Sin  n  •  W 

sin  I  ' 
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Expressing  equation  (9)  in  terms  of  the  independent  variable  r\  instead  of 
t  we  obtain  (see  Ref.  3) 


II 

^  ^  {R  e  sin  Y  + 
f^P 

S(  1  +  e  cos  Y  )1 

II 

^  ^  f  R  sin  Y  +  ^  [2 
P  P 

2 

cos  Y  +  e(l  +  cos  Y)] 

•  Si 

II  II 

3 

^  ^  cos  n  •  w 

pp 

3 

r  H  sin  T]  . 
jip  sin  I 

>  (11) 

d(ju 

d?i 

2 

f  -  R  cos  Y  +  S  • 
pe 

(l+~)sinYl  “  cos 
P 

.  dQ 
^  dTl 

II 

2 

r  H 

/  ^V2 

(PP) 

We  could  replace  the  equation  for  the  semi-major  axis  a  with  that  for  the 
semi-latus  rectum  p: 


dp  Zr^K  c 

d?f  =  p 


(12) 


V.  SMALL  ECCENTRICITY  AND  INCLINATION 

Equations  (11)  are  singular  when  e  =  0  .  Therefore  as  in  Ref.  2, 
p.  287,  we  introduce  the  variables 


H  =  e  sin  fi) 
K  =  e  cos  uo 


8 


and  replace  the  equations  in  (11)  for  e  and  u)  by 


dl^ 

dfr 


dK 

d?\ 


2  _  2 
^  f—  R  cos  T1  +  S  —  H(l  +  cos’i') 

M-  P 

+  2S^sinrf  +  S(l— ^)sinYcos'tpl+  K(l—  cosi) 

2  2 
[r  sin  rf  +  S  —  K(  1  +  cos  Y) 

M*  P 

+  2S^cosri  —  S(l— ^)sinY  sin^uol  ~  H(l— cosI) 


y  (14) 


We  used  the  variable  instead  of  ir  in  (13)  to  ameliorate  the  singu¬ 
larity  that  occurs  when  I  =  0^  or  180^  ,  which  is  also  why  we  used  rf 
instead  of  r\  as  the  independent  variable.  If  I  is  not  near  0^  or  180^  , 
say  10^  <  I<  170^  ,  we  can  employ  the  equations  in  (11)  for  I  and  Q. 
However,  if  I  is  near  0^  or  180^  ,  we  follow  Ref.  2,  p.  288  and  make  the 
change  of  variables. 


P  =  tan  I  sin  Q 
Q  =  tan  I  cos  Q 


and  replace  the  equations  in  (11)  for  I  and  Q  by 


dP 

3 

r  K 

pp 

dQ 

r\ 

drt 

pp 

W 


r  ^  4  (cos  I  —  1)  cos  Q  sin  r\  1 

^  cos  I  -* 


W 


os  ^  +  (1  ~  cos  I)  sin  0  sin  r\ 


2, 

cos  I 


•] 


(16) 


We  cannot  entirely  replace  I  ,  Q  by  P  ,  Q  because  P  ,  Q  and  the  equa¬ 
tions  they  satisfy  are  singular  when  I  =  90^  . 
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VI.  ORBIT  AVERAGED  VARIATION 


Changing  notation,  let 


=  a 
=  H 
8^  =  K 

=  !’ 


s, 


if 

I  near 

0°  or 

180 

otherwise 

if 

I  near 

0°  or 

180 

otherwise 

(17) 


Then  equations  (11),  (14)  and  (16)  can  be  written  in  the  symbolic  form 

.k 


dB 

drf 

dt 

d^ 


rl^/nl  ~ 

f  (6  . 6  ,  r|,  t) 


f^(B^ . 8^,  rf,  t) 


k  =  1,....  5 


(18) 


Exact  integration  of  these  equations  would  give  the  exact  motion  of  the  body. 

The  osculating  orbital  elements  are  slowly  varying  functions  of  time, 
unlike  the  cartesian  coordinates.  As  a  function  of  ,  the  behavior  of  an 
osculating  element  might  be  as  in  Fig.  1.  During  one  revolution  (rf  increas¬ 
ing  by  2tt)  the  value  of  an  element  oscillates;  but  from  one  revolution  to  the 
next  there  is  a  secular  or  long  periodic  trend.  If  it  is  sufficient  to  know  this 
mean  behavior  of  the  orbital  elements  rather  than  their  exact  behavior,  we 
can  follow  Ref.  3  and  use  a  method  of  orbit  averaging  which  results  in  con¬ 
siderable  savings  in  computer  time. 

1  5 

Namely,  in  a  first  order  theory  we  hold  the  elements  (B  ,  .  .  .  ,  6  ) 
constant  in  the  right  side  of  (18)  and  evaluate  t  from  rf  using  these  ele¬ 
ments  and  elliptic  orbit  formulas.  Then  the  first  order  mean  changes  in  the 
elements  during  revolution  j  are 
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Fig.  1.  Possible  behavior  of  an  osculating  elliptic  orbital  element 
as  a  function  of  the  angular  variable  rj'. 
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AB 


period 


^  =  f 

J  Ml 

_  r 

"  y?Ti 


.  .  .  ,  ri,  t)  dri  ,  k  =  1 . 5 


f^(B\  .  . .  ,  B'’,  rf,  t)  dri 


(19) 


where  rf,  is  the  value  of  rf  for  which  an  orbit  is  considered  to  commence, 
1  ^  k  k  k 

where  r)^  =  +  2it  and  when  B  ~  1  ’  value  of  B  at  the  end 

of  the  previous  revolution  (j-1)  . 


We  can  evaluate  the  integrals  in  (19)  to  second  order  by  taking 
,k 


B 


-  B^ 

-  Sj.i 


rsj  V 

(ri-Ti,)  , 

Ztt  j-1 


k  =  1,  .  .  .  .  5 


(20) 


where  ABj  ^  is  the  variation  in  the  elements  extrapolated  from  one  or  more 
immediately  preceding  revolutions.  The  first  revolution  for  which  there  is 
no  previous  revolution  can  be  integrated  twice,  the  first  time  with  the  orbital 
elements  constant  and  the  second  time  with  the  variation  in  the  orbital  ele¬ 
ments  calculated  in  the  first  integration.  Thereafter,  each  revolution  need 
only  be  integrated  once  with  the  variation  in  the  elements  in  the  integrands 
being  extrapolated  from  the  variation  determined  in  one  or  more  previous 
revolutions . 

Instead  of  using  the  orbit  averaged  first  order  linear  variation  of  the 
elements  in  the  second  order  integration,  better  accuracy  could  be  obtained 
by  using  non-orbit  averaged  first  order  variations  in  the  elements.  These 
would  have  to  be  derived  analytically. 

It  is  desirable  to  use  a  second  order  evaluation  of  the  integrals  (19) 
in  order  to  calculate  the  period  accurately  from  one  orbit  to  the  next.  How¬ 
ever,  the  accurate  calculation  of  the  secular  and  long  period  behavior  of  the 
1  5 

elements  (p  ,  .  .  .  ,  0  )  giving  the  shape  and  orientation  of  the  orbit  are  not 
so  dependent  on  the  use  of  a  second  order  theory,  except  insofar  as  their  pre¬ 
diction  might  get  out  of  phase  with  time. 

Our  main  concern  is  with  earth  satellite  orbits  for  which  the  moon 
moves  appreciably  during  an  orbital  revolution  and  the  ratio  of  the  dis¬ 
tances  of  the  satellite  and  the  moon  from  the  earth  is  not  small.  These 
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facts  cause  difficulties  with  analytic  treatments.  Therefore,  the  doctrine 
expounded  in  this  note  is  to  evaluate  the  integrals  in  (19)  numerically.  An 
accurate  numerical  second  order  orbit  averaged  theory  might  still  require 
an  analytic  non- orbit  averaged  first  order  theory  instead  of  the  first  order 
theory  given  by  (20).  It  would  be  easy  to  add  such  analytic  expressions  to 
a  computer  program  which  numerically  evaluates  the  integrals  (19).  We 
therefore  do  not  persue  the  possibilities  of  analytic  non-orbit  averaged  first 
order  theories  in  this  note. 

A  computer  program  would  run  faster  with  completely  closed  analytic 
formulas  than  it  would  if  the  integrals  in  (19)  were  evaluated  numerically. 
However,  the  numerical  integration  of  the  definite  integrals  in  (19)  is  still 
faster  than  numerically  integrating  exact  equations  of  motion  (1).  Since  the 
cartesian  coordinates  enter  implicitely  in  the  right  side  of  (1)  a  predictor- 
corrector  numerical  integration  method  is  employed  and  from  100  to  200 
or  more  steps  made  per  orbital  revolution.  The  number  of  steps  cannot  be 
decreased  even  if  less  accuracy  is  required  because  the  integration  becomes 
unstable.  On  the  other  hand,  everything  in  the  integrands  in  (19)  are  ex- 
plicite  functions  of  the  independent  variable.  Therefore,  an  integration  tech¬ 
nique  such  as  Gaussian  quadrature  can  be  employed  with  many  fewer  steps 
per  orbital  revolution. 

VII.  ELLIPTIC  ORBIT  FORMULAS 

15/^ 

Given  values  of  (B  ,  .  .  .  ,  B  )  and  rf  we  must  evaluate  the  various 
quantities  which  enter  in  the  numerical  integration  of  the  right  sides  of  (19). 
Of  course,  B^  =  a  .  By  (13)  and  (17)  we  have 


(21) 


e 


If  e  =  0  we  take  u)  =  0  and  u)  =  Q  by  convention.  If  e  >  0  we  have 


cos  u) 


sin  0) 


(22) 
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We  either  have  I  = 


.  n  =  or  by  (15) 


tan  I 
sin  Q 
cos  n 


r  4  2 

[(B*)  +  (B’)  J 


=  B  /tan  I 
=  B^/tanl 


(23) 


We  then  calculate 


Y  =  ^  -  w 

Ti  =  rf  —  0 

p  =  a(l  -  e^) 

1/2  -3/2 
n  =  fJL  '  a  ' 


(24) 


1  +  e  cos  Y 

We  calculate  the  transformation  matrix 


^11 

= 

COS  Q  cos  T]  - 

^12 

= 

—  cos  n  sin  r\ 

^13 

sin  Q  sin  I 

^21 

= 

sin  n  cos  T)  + 

^22 

= 

—  sin  Q  sin  r\ 

^23 

= 

—  cos  n  sin  [ 

^31 

= 

sin  r\  sin  I 

®32 

cos  T]  sin  I 

^33 

= 

cos  I 

(25) 


Then  the  relations  between  the  unit  vectors  e 


r 


G2»  in  the  x 


2  3 

X  ,  X 


directions  and  the  unit  vectors  e  ,  s  ,  e  defined  in  Section  III  are 

r  s  w 
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G 

r 


G 

S 


G 


w 


E  B., 
j=i  j 


Z  B.-  e. 
j=l  J 


J 


Z  B.,  s. 
i=l  J 


(26) 


and  the  cartesian  coordinates  of  the  satellite  relative  to  the  earth  are 


X  =  rB 


kl 


k  =  1,  2,  3 


(27) 


To  calculate  the  time  t  from  the  epoch  when  rf  was  ,  we  perform 
the  following  steps; 


§  =  2  tan 


5  =2  tan 


-  e  .  Y  ' 

1~  2 


M  =  §  —  e  sin  5 

o  o 

M  =  §  —  esin§ 


(28) 


t  = 


-  (M  -  M  ) 
n  o 


The  next  two  sections  derive  the  expressions  for  the  acceleration  on 
the  satellite 


1  2  3 

F  =  F  +  F  +  F 


(29) 


due  to  the  sun,  moon  and  earth  gravitational  potential  harmonics  at  a  given 


1  5 


time.  The  components  of  F  in  the 


directions  are  then 


s 


and  e 


w 


R 


F 


e 


r 


S 


F 


e 


s 


(30) 


W 


F 


e 


w 


We  define 


k 

X 

me 


k 

X 

cs 


k 

X 

se 


M 

s 


M 

m 


M 

e 


M 

c 


y 


We  have 


k 

X 

s  e 


f  DUE  TO  THE  SUN  AND  MOON 

th 

k  coordinate  of  the  moon  relative  to 
the  earth  (k  =  1,  Z,  3) 

th 

k  coordinate  of  the  earth-moon  barycenter 
relative  to  the  sun  (k  =  1,  2,  3) 

th 

k  coordinate  of  the  sun  relative  to  the 
earth  (k  =  1,  2,  3) 

mass  of  sun 

mass  of  moon 

mass  of  earth 

M  +  M 
e  m 

gravitational  constant  (so  that  \i  =  yM^) 


We  define 


k 

X 

me 


k 

X 

cs 


k 

X 

se 


M 

s 


M 

m 


M 

e 


M 

c 


y 


We  have 


k 

X 

s  e 


c 
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Given  the  time  we  read  a  magnetic  tape  and  interpolate  to  determine 
k  k 

(k  =  1,  Z,  3)  .  The  already  existing  magnetic  tape  and  computer 

subroutines  which  we  possess  give  the  coordinates  of  the  moon  and  sun  in  the 
coordinate  system  referred  to  the  mean  equinox  and  equator  of  1950.0  .  We 
have  so  far  left  the  coordinate  system  in  which  we  integrate  the  equations  of 
motion  unspecified.  Given  the  subroutines  which  already  exist,  we  shall  use 
the  inertial  system  referred  to  the  mean  equinox  and  equator  of  1950.  0  . 

Having  determined  the  position  of  the  sun  and  moon  relative  to  the 
earth,  the  perturbing  accelerations  due  to  these  bodies  on  the  motion  of  the 
satellite  relative  to  the  earth  are  the  differences  in  the  accelerations  on  the 
satellite  and  earth: 


F 


k 


V  r 
a=m,  s 


M 


a 


r  k 
ab 


ab 


k  1 

X 

ae 

3 

r 

ae 


k  =  1,  2,  3  (‘32) 


where  for  a  =  m,  s  and  A.  =  b,  e 


k 


X 


ab 


r 


aA. 


1,  2,  3 


IX.  ACCELERATION  DUE  TO  EARTH  GRAVITATIONAL  POTENTIAL 
HARMONICS 

12  3 

Let  (y  ,  y  ,  y  )  be  a  coordinate  system  with  origin  at  the  center  of 

3  . 

mass  of  the  earth,  with  y  axis  pointing  along  the  axis  of  rotation  of  the  earth, 

13  3 

with  y  axis  normal  to  y  lying  in  the  plane  which  contains  the  y  axis  and 

2 

Greenwich,  and  with  y  axis  completing  the  right  hand  system.  The  relation 

12  3 

between  this  coordinate  system  and  the  one  (x  ,  x  ,  x  )  referred  to  the  mean 
equinox  and  equator  of  1950.  0  is 


17 


k 

y  = 


k 

X 


k  =  1,  2,  3 


(33) 


where  the  orthogonal  rotation-nutation-precession  matrix  (A..)  is  discussed 
in  Section  X. 

We  introduce  polar  coordinates  (r,  0,  0)  rotating  with  the  earth  by 


We  have 


y  =  r  cos  9  cos  0 


0  ^  r  < 


y  =  r  sin  0  cos  0 
y  =  r  sin  0 


O^0<2-n- 
2  ®  2 


3  3 

^  ,  1,2  ^  ,  1,2 

Z  (x  )  =  E  (y  ) 


l=l 


1  =  1 


sin  0 


-  E  A  x"^ 
^  1  =  1 


COS  0 


/t 

V  1  —  sin  0 


:os  0 


sin  0 


— ^  2  A,, 

r  cos  0  ^  It 


_ l _  E  A  x”^ 

rcos0  Zl, 


(34) 


(35) 


(36) 


(37) 
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Outside  the  earth  the  gravitational  potential  of  the  earth  can  be  expressed 
in  spherical  harmonics  by 


U 


-  -  { 

00 

1  ~  E 

r  ( 

n=2 

00 

^  1 
y:  — 

+  E 

n=2 

h  =  l  r 

•  ^nh 

(sin  0)  > 

J 

—  P  (sin  0) 
n  n 


cos  h  9 


+  sin  h  6 


(38) 


where  and  P  ^  are  the  Legendre  ploynomials  and  generalized  Legendre 

functions,  respectively  (see  Section  XI).  The  summation  starts  with  n  =  2 
rather  than  with  n  =  1  because  the  origin  of  the  coordinate  system  is  at  the 

center  of  mass  of  the  earth. 

3 

Since  to  a  high  degree  of  approximation  the  y  axis  is  a  principal  mo¬ 
ment  of  inertia  axis,  we  have 

C^i  =  0  ,  =  0  (39) 


The  J  are  called  the  zonal  harmonic  coefficients  and  the  C  ,  ,  S  ,  are 

n  nn  nil 

called  the  tesseral  harmonic  cosine  and  sine  coefficients.  C  and  S  are 

nn  nn 

also  known  as  sectorial  harmonic  coefficients.  The  sign  and  notation  con¬ 
ventions  are  those  adopted  by  the  Smithsonian  Astrophysical  Observatory; 
see  Ref.  4. 

The  cosines  and  sines  of  multiples  of  the  longitude  9  can  be  calculated 


by 


cos  2  9 
sin  2  9 


cos  h  9 
sin  h  9 


cos  9  —  sin  9 
2  sin  9  cos  9 


>  (40) 


cos(h  —  1)  9  cos  9  —  sin(h—  1)  9  sin  9 
sin(h  —  1)9  cos9  +  cos (h—  1)9  sin0 
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The  acceleration  on  the  satellite  due  to  the  gravitational  attraction  of 

th 

the  earth  is  —  grad  U  .  Thus,  the  k  component  of  the  acceleration  in  the 
1  2  3 

(x  ,  X  ,  X  )  coordinate  system  is 


|JIX 

~T 

r 


.  ^  ‘^n  r(n+l)x^  r:i  /  •  9sin(^1 

^  -i^L  r  >■  — irj 

n=Z  r  9x 


00  n 


1 


+  |j.  Z  E  ^2 

n=2  h=l 


Cnh  =os  h  6  +  sin  h  0 


)x' 


P_^h(*in0) 


+  P^^fsiniS)  +  >>[- 


Pnh<=“»l  •  'TT 

ax 


(41) 


where  by  (36) 


a  sin  0  ^  X  .  ^ 

r  =  A3j^  -  —  sm  0 

ax 


(42) 


and  by  (37) 


—  sin  0 


ae 

ax 


A,,  - 


1  6  cos  0 


r  cos  0  "^Ik  cos  0  k 

a  X 


cos  0  —  — cos  6 


cos  9  — 

ax 


A.,  - 


1  6  cos  0 


r  cos  0  "^Zk  cos  0  k 

O  X 


sin  9  — 


k 

sin  0 
r 


Multiplying  the  first  equation  by  —  sin  9  and  the  second  by  cos  9  and  adding 
we  obtain 


S  0 


Bx 


IS  0  [■' 


k  -  |A2k‘=°=  «  -  Aik  ' 


in  sj 


(43) 
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k  3 

For  use  in  (1)  and  (30)  the  —  fix  / r  term  should  be  dropped  from 
(41).  The  order  of  the  harmonics  to  be  employed  depends  on  the  case  being 
considered.  For  all  satellites  it  is  necessary  to  include  .  For  most 
satellites  the  inclusion  of  a  few  additional  zonals  would  be  sufficient.  For  a 
synchronous  satellite  we  would  certainly  want  to  include  C_~  and  S_~  . 


X.  ROTATION  -  NUTATION  -  PRECESSION  OF  THE  EARTH 

The  transformation  matrix  (A..)  appearing  in  (33)  is  defined  by 


■^1 

^12 

■^13  ' 

■^21 

■^22 

^23 

= 

SNP 

(44) 

^31 

^32 

^3  3 

where  S  gives  the  transformation  between  coordinates  (y  .  y  ,  y  )  rotating 
with  the  earth  and  those  referred  to  the  true  equinox  and  equator  of  date, 
where  N  gives  the  transformation  between  coordinates  referred  to  the  true 
equinox  and  equator  of  date  and  those  referred  to  the  mean  equinox  and  equa¬ 
tor  of  date,  and  where  P  gives  the  transformation  between  coordinates  refer- 

12  3 

red  to  the  mean  equinox  and  equator  of  date  and  those  (x  ,  x  ,  x  )  referred 
to  the  mean  equinox  and  equator  of  1950.0  .  Matrix  multiplication  follows 
the  usual  row  x  column  rule,  so  that,  for  instance, 


(NP).. 

U 


3 

E 

k=l 


N. 


ik 


(45) 


Following  Ref.  5,  p.  482,  we  define  the  angles 


z 


9 


2304:’948T 

2304:'948T 

2004;'255T 


+  0;'302T^ 
+  l'.'093T^ 

-  01'426T^ 


+  0:'0179T^ 
+  01'0192T^ 
-  01'0416T^ 


(46) 
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where  T  is  measured  in  tropical  centuries  of  36524.21988  ephemeris  days 
from  the  epoch  1950.0  (JED  2433282.423)  to  the  instant  of  interest.  Then 
by  Ref.  6,  p.  3  1,  the  precession  matrix  P  at  this  instant  is 


^11 

COS  C 
'  0 

cos  9  cos  z  — 

sin  C  sin  z 
^0 

^12 

= 

—  sin  C 

cos  9  cos  z 
0 

—  cos  C  sin  z 
^0 

^13 

= 

—  sin  9 

cos  Z 

^21 

= 

cos  C 
^0 

cos  9  sin  z  + 

sin  C  cos  z 
^0 

^22 

—  sin  C 

cos  9  sin  z 
0 

+  cos  C  cos  z 
^0 

^23 

—  sin  9 

sin  z 

P3I 

= 

cos  C 
^0 

sin  9 

^32 

= 

-  sin  Q 

sin  9 

0 

^33 

cos  9 

Let  T  denote  the  time  from  the  epoch  1950.0  (JED  2433282.423)  in  units 
of  10,000  ephemeris  days.  Then  by  Taylor’s  theorem  we  have 


E 

n  =  0 


J_ 

nl 


^  P-1, 
Jk 

dT 


n 


T 


T  =  0 


j,  k  =  1,  2,  3 

(48) 


Treating  the  coefficients  in  (46)  as  exact,  some  simple  calculations  show 
that  the  terms  up  to  the  fifth  power  in  the  Taylor  expansions  (48)  are: 


22 


12 


13 


21 


22 


23 


31 


32 


1.0  -  2. 22  603  398052517  x  lO'^T^  -  2.  6903  5  853253  66  x  lO"*^!^ 
+  8.  191221606878  x  IO’^^t'^  +  1.  79948222850  x  10‘^^t^ 

-  6.  119064710033514  x  lO'^T  -  5. 06975739290688  x  10‘^t^ 

+  4.  5321716219079  x  lO'^T^  +  8. 619581795926  x  10‘^^t^ 

-  1.02943658327  x  10‘^^t^ 

-2.660399722772102  x  10‘^t  +  1.54818397804898  x  lO'^T^ 

+  1.  9729201591810  X  10"®t^  +  1.  960730253191  x  IO'^^t"^ 

-  4.  392983  54075  x  lO'^'^T^ 

6.  119064710033514  x  10'^t  +  5.06975739290688  x  10’'^t^ 

-  4.  5321716219079  x  10"®t^  -  9.  636891635856  x  IO'^^t'^ 

+  1.02604298897  x  10‘^^t^ 

1.0  -  1.  87214764627888  x  10‘^t^  -  3.  1022173551368  x  10'^t^ 
+  6.  88247882  553  5  x  10‘^  +  1.  91215207447  x  lO'^'^T^ 

-  8. 13957902909886  x  10‘^t^  -  5. 8309700675934  x  IO'^^t^ 

+  2.  994360606802  x  IO’^^t'^  +  5.  71739459043  x  IO'^^t^ 

2.  660399722772102  x  lO'^T  -  1. 54818397804898  x  10‘^t^ 

-  1. 9729201591810  x  10‘®t^  +  3. 791379581151  x  IO’^^t^ 

+  4.  50404085077  x  lO’^'^T^ 


(49) 


A  ?  1 D 

-  8.  13957902909886  x  10  t  +  1.81682  68497009  x  lO'  t 

+  3.  024323052660  x  10‘^^t'^  +  2.  58550054981  x  IO'^^t^ 

1.0  -  3.  53886334246294  x  10‘^t^  +  4.  1187882260017  x  IO'^^t^ 
+  1.308742781343  x  IO'^^t'^-  1.  12669845971  x  10‘^^t^ 


The  mean  obliquity  of  the  ecliptic  is  (see  Ref.  6,  p.  98) 


e  =  23°27'08:’26  -  46:'845T  -  0:'0059T^  +  0:'00181T^ 


(50) 
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where  T  is  measured  in  Julian  centuries  of  36525  ephemeris  days  from 
the  epoch  1900  January  0.  5  E.  T.  =  JED  2415020.  0  to  the  instant  of 
interest.  Let  AY  and  Ac  be  the  nutations  in  longitude  and  obliquity,  re¬ 
spectively,  as  given  by  the  series  in  Ref.  6,  p.  44-45.  The  true  obliquity 
of  the  ecliptic  is  then 


e  =  e  +  Ac 
o 


(51) 


The  nutation  matrix  is  given  by  (see  Ref.  6,  p.  43) 


’^11 

^12 

1 - 

CO 

1 

—  ay  cos  e 

—  AY  sin  e 

N  = 

^21 

^22 

^23 

= 

ay  cos  e 

1 

-  Ae 

1 - 

"^32 

^33. 

AY  cos  e 

Ae 

1 

Sidereal  time  is  defined  as  the  hour  angle  of  the  first  point  of  Aires 
(equinox  point).  We  therefore  have 


‘^11 

^12 

^13 

COS  0 

sin  © 

0 

S  = 

^21 

^22 

^23 

= 

—  sin@ 

cos  © 

0 

_S31 

^32 

^33_ 

0 

0 

1 

where  0  is  the  Greenwich  true  sidereal  time.  By  Ref.  6,  p.  72, 

0  =  0^  +  AY  cos  e  (54) 

where  0  is  the  Greenwich  mean  sidereal  time.  By  Ref.  5,  p.  478,  the 
o  _ 

Greenwich  mean  sidereal  time  0  at  0  universal  time  on  the  day  of  inter- 

o 

est  is 
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(55) 


Q  =  6^38^45^836  +  8,  640,  184®542T  +  o!o929T^ 
o 


where  T  denotes  the  number  of  Julian  centuries  of  3  6525  days  which,  at 
midnight  beginning  of  day,  have  elapsed  since  mean  noon  on  1900  January  0 
at  the  Greenwich  meridian  (Julian  Date  2415020.0).  The  Greenwich  mean 
sidereal  time  0  ^  at  a  given  instant  UT  of  universal  time  on  that  day  is  then 


©0=00+ 


d© 

c 

dt 


UT 


(56) 


where  by  Ref.  6,  pp.  75-76, 
d@ 

=  (1.002737909265  +  0.  589  x  10"  ^T) 

sidereal  time  seconds  per  universal  time  second  (57) 

The  equations  of  motion  are  integrated  as  a  function  of  Ephemeris  Time  ET, 

and  the  value  of  ET  —  UT  is  given  in  Ref.  5,  p.  vii  . 

It  was  implicitely  assumed  in  Section  IX  that  the  coordinate  system 
12  3. 

(y  >  y  »  y  )  defined  in  terms  of  the  axis  of  rotation  of  the  earth  is  fixed  in 
the  crust  of  the  earth.  This  is  not  exactly  true  and  in  fact  if  (44)  is  to 
define  the  transformation  matrix  for  coordinates  fixed  in  the  earth  to  those 
referred  to  the  mean  equinox  and  equator  of  1950.0  it  should  read  A  =  WSNP  . 
However,  the  wobble  matrix  W  is  so  nearly  the  identity  matrix  that  it  can  be 
ignored  for  the  problem  discussed  in  this  note. 

XI.  LEGENDRE  POLYNOMIALS  AND  FUNCTIONS 

By  Ref.  7,  pp.  83  and  327,  the  definitions  of  the  Legendre  polynomials 
P^  and  generalized  Legendre  functions  P^j^  are 
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(58) 


Po(Z) 

P^(Z) 

P  (Z)  =  P^(Z) 

no  n 

p„h(z)  =  P„|Z) 


2”n! 


d”(Z^  -  1)^ 


dZ 


n 


=  1,  2, 


n  =  0,  1 ,  2, .  .  . 

h  =  1 . n 


(59) 


From  these  definitions  it  follows  that 


Pn(Z) 


[?] 

E  (-1)^ 


(2n-2i)  ! 


,n-2i 


i=0 


2  (n-i)  !  (n-2i)  !  i  1 


(60) 


r— 1 


(2n-2i)  ! 


,n-h-2i 


.  n 


(61) 


i  =  0  2  (n-i)  1  (n-h-2i)  I  i  I 


where  [x]  denotes  the  largest  integer  less  than  or  equal  to  x  .  A  compu¬ 
ter  subroutine  to  evaluate  the  Legendre  polynomials  and  functions  and  their 
derivatives  should  use  the  various  recursion  formulas  that  exist. 

In  Ref.  4  expansion  (38)  is  employed  with  and 

in  place  of  C  ,  ,  S  ,  and  P  ,  ,  where 
nn  nn  nn 


p  =  ^  •  p  (6: 

nh  \  (n  +  n)  [  nn 

This  is  because  the  integral  of  P^j^(sin  0)  times  cos  h  6  or  sin  h  0  over 
the  sphere  is  4ir  .  However,  Ref.  4  does  not  use  the  corresponding  normali¬ 
zation  for  the  zonal  harmonics 


P  =  y'Zn  +  1  P 
n  n 
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(63) 


In  using  the  values  of  given  in  Ref.  4  from  fits  to  satellite 

data  these  different  definitions  should  be  taken  into  account. 


XII.  GAUSSIAN  QUADRATURE 

The  method  of  numerical  integration  that  should  be  used  in  evaluating 

th 

the  definite  integrals  (19)  is  Gaussian  quadrature.  Let  be  the  n 

Legendre  polynomial.  Following  Ref.  8,  pp.  319-325,  we  have 


f(y)  dy 


m 

S  H  f(y  ) 
i  =  l  ^  ^ 


where  y.  is  the 


.th 

1 


zero  of 


P  (y) 

m 


and  where 


(64) 


H. 

1 


m  P  (y  )  P;^(y  ) 

m- 11  mi 


2(1  -  y^) 


The  formula  (64)  is  exact  for  polynomials  f  of  order  2m—  1  .  For  the 
integrals  of  interest  to  us  we  have 


/  f(Tl)  dTt  =  S  I  ^ 
•'tT,  j  =  l  J2j(] 


1 


f(?T)  diT 


(ill)  + 


since  rf^  “  •  To  evaluate 


I 


m  dTi 
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we  define 


/b  a^  1  ^b'4*a^ 

T)  =  (  2  )  y  +  (  2  ^ 

d?r  =  (•^^)  dy 

so  that 

f. 

* 

f(Tf)  d?|  =  f  [(^^)  y  +  (^-|-^)]  dy 

b  — a  ^  TTrf/b— a,  .  ,b  +  a.l 

=  2  2  '  n  +  <  2  >J 

1=1  L  -» 

Tables  of  ,  H^(i  =  1,  .  .  .  ,  m)  are  given  in  Ref.  9  for  m  =  1,  .  .  .  ,  16 
The  values  of  t  and  m  should  be  input  and  chosen  to  give  the  desired  accu¬ 
racy  in  the  fastest  time. 
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